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Abstract 



A general, iterative, method for the description of evolving self-gravitating 
relativistic spheres is presented. Modeling is achieved by the introduction 
of an ansatz, whose rationale becomes intelligible and finds full justification 
within the context of a suitable definition of the post-quasistatic approxima- 
tion. As examples of application of the method we discuss three models, in 
the adiabatic case. 
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I. INTRODUCTION 



The problem of general relativistic gravitational collapse has attracted the attention of 
researchers since the seminal paper by Oppenheimer and Snyder |L]. The motivation for 
such interest is easily understood: the gravitational collapse of massive stars represents one 
of the few observable phenomena, where general relativity is expected to play a relevant 
role. Ever since that work, much was written by researchers trying to provide models of 
evolving gravitating spheres. However this endeavour proved to be difficult and uncertain. 
Different kinds of advantages and obstacles appear, depending on the approach adopted for 
the modeling. 

Thus, numerical methods (see || and references therein) are enabling researchers to 
investigate systems which are extremely difficult to handle analytically. In the case of 
General Relativity, numerical models have proved valuable for investigations of strong field 
scenarios and have been crucial to reveal unexpected phenomena ||. Even specific difficulties 
associated with numerical solutions of partial differential equations in presence of shocks are 
being overpassed Q. By these days what seems to be the main limitation for numerical 
relativity is the computational demands for 3D evolution, prohibitive in some cases ||. 
Nevertheless, purely numerical solutions usually hinder to catch general, qualitative, aspects 
of the process. 

On the other hand, analytical solutions although more suitable for a general discussion 
(see and references therein), are found, either for too simplistic equations of state and/or 
under additional heuristic assumptions whose justification is usually uncertain. 

Therefore it seems useful to consider nonstatic models which are relatively simple to 
analyze but still contain some of the essential features of a realistic situation. 

Accordingly, it is our purpose in this work to present an approach for modeling the 
evolution of self-gravitating spheres, which may be regarded as a "compromise" between 
the two approaches mentioned above (analytical and numerical). 

Indeed, the proposed method, starting from any interior (analytical) static spherically 
symmetric ( "seed" ) solution to Einstein equations, leads to a system of ordinary differential 
equations for quantities evaluated at the boundary surface of the fluid distribution, whose 
solution (numerical), allows for modeling the dynamics of self-gravitating spheres, whose 
static limit is the original "seed" solution. 

The approach is based on the introduction of a set of conveniently defined "effective" 
variables, which are effective pressure and energy density, and an heuristic ansatzs on the 
latter [0], whose rationale and justification become intelligible within the context of the 
post-quasistatic appproximation defined below. In the quasistatic approximation (see next 
section), the effective variables coincide with the corresponding physical variables (pressure 
and density) and therefore the method may be regarded as an iterative method with each 
consecutive step corresponding to a stronger departure from equilibrium. In this work we 
shall restrain ourselves to the post-quasistatic level (see next section for details). 

At this point it is important to stress a crucial difference between this method and 
the one proposed many years ago with a similar structure, but based on radiative Bondi 
coordinates (see || and references therein): in the latter the introduced effective variables 
do not coincide with the corresponding physical variables in the quasistatic approximation 
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(they do coincide in the static limit), and accordingly the ansatz on those variables remains 
as an heuristic assumption, only justified by the eventual suitability of the obtained models. 

The fluid distribution under consideration will be assumed to be dissipative. Indeed, 
dissipation due to the emission of massless particles (photons and/or neutrinos) is a char- 
acteristic process in the evolution of massive stars. In fact, it seems that the only plausible 
mechanism to carry away the bulk of the binding energy of the collapsing star, leading to a 
neutron star or black hole is neutrino emission ||. Consequently, in this paper, the matter 
distribution forming the self gravitating object will be described as a dissipative fluid. 

In the diffusion approximation, it is assumed that the energy flux of radiation (as that of 
thermal conduction) is proportional to the gradient of temperature. This assumption is in 
general very sensible, since the mean free path of particles responsibles for the propagation 
of energy in stellar interiors is in general very small as compared with the typical length of 
the object. Thus, for a main sequence star as the sun, the mean free path of photons at the 
centre, is of the order of 2 cm. Also, the mean free path of trapped neutrinos in compact 
cores of densities about 10 12 g.cm.~ 3 becomes smaller than the size of the stellar core []10, 11 



Furthermore, the observational data collected from supernovae 1987A indicates that the 
regime of radiation transport prevailing during the emission process, is closer to the diffusion 
approximation than to the streaming out limit W2 . 



However in many other circumstances, the mean free path of particles transporting en- 
ergy may be large enough as to justify the free streaming approximation. Therefore our 
formalism will include simultaneously both limiting cases of radiative transport (diffusion 
and streaming out), allowing for describing a wide range situations. 

Besides the usual physical variables (energy density, pressure, velocity, heat flow, etc.) 
we shall also incorporate into discussion other quantities which are expected to play an 
important role in the evolution of evolving self-gravitating systems, such as the Weyl tensor, 
the shear of the fluid and the Tolman mass. Therefore these quantities will be calculated 
and used in the process of modeling. It is also worth mentioning that although the most 



common method of solving Einstein's equations is to use commoving coordinates (e.g. |13 
||), we shall use noncomoving coordinates, which implies that the velocity of any fluid 
element (defined with respect to a conveniently chosen set of observers) has to be considered 
as a relevant physical variable ( Hl4|). 

The plan of the paper is as follows. In Section 2 we define the conventions and give the 
field equations and expressions for the kinematical and physical variables we shall use, in 
noncomoving coordinates. The proposed approach is presented and explained in Section 3. 
In Section 4 we illustrate the method by means of three examples. Finally a discussion of 
results is presented in Section 5. 

II. RELEVANT EQUATIONS AND CONVENTIONS 

A. The field equations 

We consider spherically symmetric distributions of collapsing fluid, which for sake of 
completeness we assume to be locally anisotropic, undergoing dissipation in the form of heat 
flow and/or free streaming radiation, bounded by a spherical surface S. 
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The line element is given in Schwarzschild-like coordinates by 



sin 



(1) 



where v(t,r) and X(t, r) are functions of their arguments. We number the coordinates: 
x° = t; x 1 = r; x 2 = 8; x 3 = 0. 

The metric ([I]) has to satisfy Einstein field equations 



which in our case read UTS 
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where dots and primes stand for partial differentiation with respect to t and r, respectively 
In order to give physical significance to the components we apply the Bondi approach 
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Thus, following Bondi, let us introduce purely locally Minkowski coordinates (r, x, y, z) 

dr = e u / 2 dt ; dx = e x ^ 2 dr ; dy = rdd ; dz = rsinOdcf). 
Then, denoting the Minkowski components of the energy tensor by a bar, we have 



T 1 = T 1 ■ % 
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Next, we suppose that when viewed by an observer moving relative to these coordinates 
with proper velocity uj in the radial direction, the physical content of space consists of an 
anisotropic fluid of energy density p, radial pressure P r , tangential pressure Pj_, radial heat 
flux q and unpolarized radiation of energy density e traveling in the radial direction. Thus, 
when viewed by this moving observer the covariant tensor in Minkowski coordinates is 

/ p + e -q-e \ 
-q-e P r + e 
P ± 

\ Pi/ 
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Then a Lorentz transformation readily shows that 
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Note that the coordinate velocity in the (t,r,8,<p) system, dr/dt, is related to 
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Feeding back (fTHTOD into we get the field equations in the form 
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Observe that if v and A are fully specified, then (]14HT7|) becomes a system of algebraic 
equations for the physical variables p, P r , P±, u, Q and e. Obviuosly, in the most general 
case when all these variables are non-vanishing, the system is underdetermined, and two 
equations of state should be given. In general, whenever Q ^ a transport equation has to 
be assumed. In the case originally considered by Bondi |15| (locally isotropic fluid and free 
streaming regime, Q — 0) the system is closed. For the adiabatic (e = Q — 0), and locally 
isotropic fluid (P r = P±) the system is overdetermined, and a constraint on the physical 
variables appears. 

At the outside of the fluid distribution, the spacetime is that of Vaidya, given by 

ds 2 = ^1 - 2M ^ ^j du 2 + 2dudK - U 2 (d9 2 + sin 2 9d0 2 ) , (18) 

where u is a coordinate related to the retarded time, such that u = constant is (asymp- 
totically) a null cone open to the future and 1Z is a null coordinate (gun = 0). It should 
be remarked, however, that strictly speaking, the radiation can be considered in radial free 
streaming only at radial infinity. 

The two coordinate systems (t, r, 9, <fi) and (u, 71, 9, (p) are related at the boundary surface 
and outside it by 



u = t _ r _2Mln( K — - 1 j , (19) 

11 = r. (20) 

In order to match smoothly the two metrics above on the boundary surface r = rs(t), we 
first require the continuity of the first fundamental form across that surface. Then 



l _2M +2 dRz 
Ry, du 



du 2 , (21) 



where R = Ry,{u) is the equation of the boundary surface in (u, R, 9, 0) coordinates. 
From (gTp, using (0), flT|) and (0) it follows 



2M , . 

e us = 1 - (22) 

e- x - = 1 - — . (23) 

Where, from now on, subscript £ indicates that the quantity is evaluated at the boundary 
surface E. 

Next, the unit vector n M , normal to the boundary surface, has components 



dR 
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where + indicates that the components are evaluated from the outside of S, and (3 is given 
by 

„ / 2M(u) dR s \~ 1/2 , s 

0=1-— -Li +2-^ . (25) 



i?s du 

The unit vector normal to £, evaluated from the inside, is given by 

n { - ] = (-r S 7, 7, 0, 0) , (26) 

with 

7= (e~ As -r|e-^)" 1/2 . (27) 
Let us now define a time-like vector such that 

^+)=^ + / 3f^p^ (2g) 

and 

t^-J = — + rrf (29) 

Then, junction conditions across E, require (besides ([21]) ) 

(T^tAi")^ = {T^n v )i ] , (30) 

(T^n^") « = (T^n^p , (31) 

where the expressions for the energy momentum tensor at both sides of the boundary surface 
are 

T $ = (P + ^l) u v u v ~ P±9^u + (P r - P±) s M s„ + q^Uy + q v u^ + elj^ (32) 

and 
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where w M denotes the four velocity of the fluid, s M is a radially directed space-like vector 
orthogonal to u^, / M is a null outgoing vector, and 



q» = Q (W A ^ )/2 , 1,0, 0) 
Then it follows from ([30]) and (|3l|) 
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Eqs. (|2l|), (^1) and (|39| ) are the necessary and sufficient conditions for a smooth matching 
of the two metrics ([!]) and ([18]) on E. Combining (|38|) and (]3^) we get 
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expressing the discontinuity of the radial pressure in the presence of heat flow, which is a 
well known result [ |T6|j . 

Next, it will be useful to calculate the radial component of the conservation law 
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After tedious but simple calculations we get 
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which in the static case becomes 



2 r 



(43) 



representing the generalization of the Tolman-Oppenheimer-Volkof equation for anisotropic 
fluids flTJ. 



B. The kinematical variables 



The components of the shear tensor are defined by 
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(45) 



denote the projector onto the three space orthogonal to w M , the expansion and the four 
acceleration, respectively. 
A simple calculation gives 
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and for the shear scalar a 
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C. The Tolman mass 

The Tolman mass for a spherically symmetric distribution of matter is given by (eq.(24) 
in M) 
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where £ denotes the usual gravitational lagrangian density (eq.(10) in |L8|]). Although 
Tolman's formula was introduced as a measure of the total energy of the system, with 
no commitment to its localization, we shall define the mass within a sphere of radius r, 
completely inside S, as 



m T = 4vr f r 2 e^ +A )/ 2 (T ° - T\ - 2T 2 ) dr 
Jo 



+ 2 



2 p (v+X)/2^_ 



r e 



d£ 



dt \d[d(g^y/^)/dt] 



(56) 



This extension of the global concept of energy to a local level fl9| is suggested by the 
conspicuous role played by m<r as the "effective gravitational mass" , which will be exhibited 
below. Even though Tolman's definition is not without its problems |19|,|2(|, we shall see 
that rriT, as defined by (|56|), is a good measure of the active gravitational mass, at least for 
the systems under consideration. 

Let us now evaluate expression (|5*6"|). The first integral in that expression 



(57) 
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may be transformed to give (see ETJ for details) 
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where the mass function m, as usually is defined by 

e -A(r,t) = l _ 2m( r ,t)/r. 

Next, from (eq.(13) in pf) 
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and so the second integral (//) in (|56D may be expressed as 
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This is, formally, the same expression for in terms of m and that appears in the static 
(or quasi-static) case (eq.(25) in |E3|). 
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Replacing by (|j) and m by fl59|), one also finds 

m T = e {u - x)/2 v' —. (63) 

This last equation brings out the physical meaning of tut as the active gravitational mass. 
Indeed, it can be easily shown |23[ that the gravitational acceleration a of a test particle, 



instantaneously at rest in a static gravitational field, as measured with standard rods and 
coordinate clock is given by 

a = — = "7? • (64) 

A similar conclusion can be obtained by inspection of Eq. ( f43l) (valid only in the static 
or quasi-static case) [p^j . In fact, the first term on the right side of this equation (the 
"gravitational force" term) is a product of the "passive" gravitational mass density (p + P r ) 
and a term proportional to my/r 2 . 

D. The Weyl tensor 

Since the publication of Penrose's work there has been an increasing interest in 
studying the possible role of Weyl tensor (or some function of it) in the evolution of self- 



gravitating systems . This interest is reinforced by the fact that for spherically symmetric 
distribution of fluid, the Weyl tensor may be defined exclusively in terms of the density 
contrast and the local anisotropy of the pressure (see below), which in turn are known to 



affect the fate of gravitational collapse p? |. 

Now, using Maple V, it is found that all non-vanishing components of the Weyl tensor 
are proportional to 



where 
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corresponds to the contribution in the static case. 

Also, from the field equations and the definition of the Weyl tensor it can be easily shown 
that (see [2T| for details) 



W = -f J\ 3 (Tffdr + fr* (Tf - T, 1 ) . (67) 



III. THE METHOD 

We have now available all the ingredients required to present our method, however before 
doing so some general considerations will be necessary. 
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A. Equilibrium and departures from equilibrium 



The simplest situation, when dealing with self-gravitating spheres, is that of equilibrium 
(static case). In our notation that means that u = e = Q = 0, all time derivatives vanishes, 
and we obtain the generalized Tolman-Oppenheimer-Volkof equation (|43|). 

Next, we have the quasistatic regime. By this we mean that the sphere changes slowly, 
on a time scale that is very long compared to the typical time in which the sphere reacts to 
a slight perturbation of hydrostatic equilibrium, this typical time scale is called hydrostatic 
time scale [^8|] (sometimes this time scale is also referred to as dynamical time scale, e.g. 
see the third reference in |]28|). Thus, in this regime the system is always very close to 
hydrostatic equilibrium and its evolution may be regarded as a sequence of static models 
linked by (|TT[). This assumption is very sensible because the hydrostatic time scale is very 
small for many phases of the life of the star. It is of the order of 27 minutes for the Sun, 4.5 
seconds for a white dwarf and 10 -4 seconds for a neutron star of one solar mass and 10 Km 
radius. It is well known that any of the stellar configurations mentioned above, generally, 
change on a time scale that is very long compared to their respective hydrostatic time scales. 
Let us now translate this assumption in conditions to uo and metric functions. 

First of all, slow contraction means that the radial velocity uo as measured by the 
Minkowski observer is always much smaller than the velocity of light (u <C 1). Therefore 
we have to neglect terms of the order 0(uj 2 ). 

Then (H3) yields 



AH = 8nre u 

2 2 



(P r + e)' + (p + P r + 2e) - - 2 



(68) 



(observe the contribution of e to, both, P r and p and the fact that e, u and Q are of the 
same order of smallness, in this approximation). 

Since by assumption, in this regime the system is always (not only at a given time t) 
in equilibrium (or very close to), (f43l) and (|68|) imply then, for an arbitrary slowly evolving 
configuration 

A « z>A « A 2 « 0, (69) 

and of course, time derivatives of any order of the left hand side of the hydrostatic equilibrium 
equation must also vanish, for otherwise the system will deviate from equilibrium. This 
condition implies, in particular, that we must demand in this regime 

v w 0. 



Finally, from the time derivative of (y), and using (|ToD , it follows that 



Co « 0(A, Acj,z>cj). (70) 

which implies that we have also to neglect terms linear in the acceleration. On purely 
physical considerations, it is obvious that the vanishing of uj is required to keep the system 
always in equilibrium. 
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Thus, in the quasistatic regime we have to assume 



0{u 2 ) = A 2 



Az> = A 



0. 



(71) 



implying that the system remains in (or very close to) equilibrium. However, during their 
evolution, self-gravitating objects may pass through phases of intense dynamical activity, 
with time scales of the order of magnitude of (or even smaller than) the hydrostatic time 
scale, and for which the quasi-static approximation is clearly not reliable (e.g., the collapse 
of very massive stars [ 29[ and the quick collapse phase preceding neutron star formation, see 
for example |3(| and references therein). In these cases it is mandatory to take into account 
terms which describe departure from equilibrium. 



B. The effective variables and the post quasistatic approximation 

Let us now define the following effective variables: 

p = To= T^ + (i-^)V 2 +e ' ( 72 ) 



In the quasistatic regime the effective variables satisfy the same equation ([43]) as the 
corresponding physical variables (taking into account the contribution of e to the "total" 
energy density and radial pressure, whenever the free streaming approximation is being 
used). Therefore in the quasistatic situation (and obviously in the static too), effective and 
physical variables share the same radial dependence. Next, feeding back ([7^) and (|73|) into 
(|i~4l) and (|15"1) , these two equations may be formally integrated , to obtain: 

m = Air r 2 pdr, (74) 



v = ^+ —-^-dr. 75 

y rs r(r - 2m) 

From where it is obvious that for a given radial dependence of the effective variables, the 
radial dependence of metric functions become completely determined. 

With this last comment in mind, we shall define the post-quasistatic regime as that 
corresponding to a system out of equilibrium (or quasiequilibrium) but whose effective vari- 
ables share the same radial dependence as the corrresponding physical variables in the state 
of equilibrium (or quasiequilibrium). Alternatively it may be said that the system in the 
post-quasistatic regime is characterized by metric functions whose radial dependence is the 
same as the metric functions corresponding to the static (quasistatic) regime. The rationale 
behind this definition is not difficult to grasp: we look for a regime which although out of 
equilibrium, represents the closest possible situation to a quasistatic evolution (see more on 
this point in the last Section). 
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C. The algorithm 

Let us now outline the approach that we propose: 

1. Take an interior solution to Einstein equations, representing a fluid distribution of 
matter in equilibrium, with a given 

Pst = p{r) P rs t = P r {r) 

2. Assume that the r dependence of P and p is the same as that of P rst and p st , respec- 
tively. 

3. Using equations ([75]) and (|74"D, with the r dependence of P and p, one gets m and v 
up to some functions of t, which will be specified below. 

4. For these functions of t one has three ordinary differential equations (hereafter referred 
to as surface equations), namely: 

(a) Equation (|i~3| ) evaluated on r = r s . 

(b) Equation (|42|) evaluated on r = r^. 

(c) The equation relating the total mass loss rate with the energy flux through the 
boundary surface. 

5. Depending on the kind of matter under consideration, the system of surface equa- 
tions described above may be closed with the additional information provided by the 
transport equation and/or the equation of state for the anisotropic pressure and/or 
additional information about some of the physical variables evaluated on the boundary 
surface (e.g. the luminosity). 

6. Once the system of surface equations is closed, it may be integrated for any particular 
initial data. 

7. Feeding back the result of integration in the expressions for m and z/, these two func- 
tions are completely determined. 

8. With the input from the point 7 above, and using field equations, together with the 
equations of state and/or transport equation, all physical variables may be found for 
any piece of matter distribution. 

D. The Surface equations 

As it should be clear from the above the crucial point in the algorithm is the system of 
surface equations. So, let us specify them now. 
Introducing the dimenssionles variables 

A = r E /m s (0), 
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F = 1 - 2M/A, 
M = m E /m E (0), 

= UJ S , 

a = t/m s (0), 

where m E (0) denotes the total initial mass, we obtain the first surface equation by evaluating 
(P"3|) at r = r s , one gets 



(76) 



Next, using junction conditions, one obtains from (]59[) , ([14]) and ([[7]) evaluating at 
r = r s , that 



da 



with 



E = 47rr|(e s + &) 



(77) 



(78) 



where the first and second term on the right of (|77|) , represent the gravitational redshift and 
the Doppler shift corrections, respectively. 

Then, defining the luminosity perceived by an observer at infinity as 

dM 



da 



we obtain the second surface equation in the form 

dF F 



da A 



1 -F)Sl + 2L/A. 



(79) 



The third surface equation may be obtained by evaluating at the boundary surface, the 
conservation law T£ = 0, which reads 



P + 



(p + P)(4:7ir 3 P + m) 
r(r — 2m) 



Anr(r — 2m) 



m + 



3m 



;,2 



r — 2m 



rrw\ 2, ~ N . 
— J+-(Pl-P). (80) 



Now, in the following section we consider two relatively simple models with a separable 
effective density, i.e., p = f(t)h(r); thus equation (p0|) evaluated at the boundary surface 
leads to 



dfl 
da 



n 2 



8F 
~A 



+ 2Fk(r s ) + 4vrp E A(3 - ft 2 



F 













E(l + Q) 
47rr| 



(81) 



where 
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R 



P 



P + p 
1 — 2m/r 



(AnrP + 



J s 



E = E(l + tt), 



and 



d 
dry 



-In 



1 



r-E, 



drr 2 h{r) / h{r^) 



(82) 



(83) 



(84) 



Before analyzing specific models, some interesting conclusions can be obtained at this 
level of generality. One of these conclusions concerns the condition of bouncing at the surface 
which, of course, is related to the occurrence of a minimum radius A. According to (76) this 
requires £1 = 0, and we have 



d 2 A 
da 2 



da 



or using 



dQ 
da 



(£1 = 0) 



F 



R 




E 



_l_£ 



47T7V 2 



(85) 



(86) 



Observe that a positive energy flux (E) tends to decrease the radius of the sphere, i.e., it 
favors the compactification of the object, which is easily understandable. The same happens 
when R > or Pj_ s < 0. The opposite effect occurs when these quantities have the opposite 
signs. Now, for a positive energy flux the sphere can only bounce at its surface when 



da 



0) > 0. 



According to 



this is equivalent to 

-R(Q : 



0) 



2P 



.4 



> 0. 



(87) 



A physical meaning can be associated to this equation as follows. For non-radiating, 
static configuration, R as defined by ([S2|) consists of two parts. The first term which together 
with — [2(Pj_ — P r )/r]s represents the hydrodynamical force (see (|43|) ) and the second which 
is of course the gravitational force. The resulting force in the sense of increasing r is precisely 
— R + [2(Pj_ — P r )/r]s, if this is positive a net outward acceleration occurs and vice-versa. 
Equation ( |37D is the natural generalization of this result for general non-static configurations. 

As mentioned before, besides the surface equations, in some cases (depending on the 
type of matter under consideration) further information has to be provided in the form of 
equation of state for the tangential stresses and/or transport equation. In the next Section 
we shall illustrate our method with three examples, one of which refers to an anisotropic 
fluid, and for which we shall further assume the equation of state (see [31], [32[0 

C(P + p)(47rr 3 P + m) 



P, ~Pr 



2m) 



(88) 



where C is a constant. 
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IV. EXAMPLES 



The only purpose of the present Section is to illustrate the proposed method. For simpli- 
fication we shall consider only the adiabatic case (e = Q — 0). For all these models we shall 
calculate the physical and geometrical variables for any piece of matter, as function of the 
timelike coordinate. In spite of the simplicity of the models, some interesting conclusions 
about the physical meaning of different variables may be reached. 

One of the models has as the "seed" solution the well known Schwarzshild interior solu- 
tion, whose properties have been extensively discussed in the literature. The second example 
is based on an anisotropic fluid without radial pressure. Models of this kind have also been 
discussed extensively since the original Einstein paper (see [0). Finally, the third example 
represents the dynamic version of the Tolman VI static solution ( [0]), whose equation of 
state, as is well known, approaches that for highly compressed Fermi gas. 



This model is inspired in the well known interior Schwarzschild solution. Accordingly we 
take 



where A; is a function of t to be determined from junction conditions ([5]), which in terms of 
effective variables, becomes 



A. Schwarzschild type model 



P = /(*). 



(89) 



where / is an arbitrary function of t 



And the expression for P is 




(90) 



(91) 




(92) 




(93) 



with 



[1 



(l-F)(r/r s ) 2 ] 1/2 



and 



X = 3(fi 2 + l)(l-F), 



V>= (3ft 2 + l)(l-F). 
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And for the metric functions m and v we get, using ( [74"D and (|75|) 



(94) 



2(1 -F) 

The third surface equation for this model becomes 



e » = J X f ^ [> . (95) 



da 2A (7-3^ + 3F(^-l)). (96) 



This equation together with ( ffq) and ( ff9|) form the set of surface equations for this 
model. We have integrated it numerically and from this integration, all physical variables 
are found for any piece of the fluid distribution, following the algorithm described above. 

Figures (l)-(5) exhibit the behaviour of p, P, uj, a and W for an initially contracting 
configuration, as function of a and different pieces of matter. 

Figure (6) shows the profile of u as function of r/rs for a = 10. 



B. Lemaitre Florides type model 

This model has as the "seed" solution a configuration with homogeneous energy density 
and vanishing radial pressure. Configurations of this kind were suggested by the first time 
by Lemaitre ( [§§])■ 

The corresponding effective variables now are: 

P = fit) (97) 

and 

P = 0. (98) 



Observe that in this case, because of (|97|) and ([98"D , it follows from fl72|) and (ff3"D that 
the radial pressure is discontinues at the boundary surface, with 

P r2 = -3(1 - F)tt 2 /8nrl, (99) 

for otherwise either ps or Q should vanish at S. Therefore the only way to "dynamize" 
this model is by relaxing boundary conditions, allowing for the presence of a kind of surface 
tension. 

Once the effective variables are defined, we only need the value of the tangential pressure 
at the boundary to close the system of surface equations. This is obtained by evaluating 
(H) at E. 

Next, following the algorithm, all physical variables may be found for any piece of ma- 
terial as functions of the timelike coordinate. Although we are not going to exhibit them 
here, because the graphics are not particularly illuminating, we wanted to present an exam- 
ple which, besides the fact that implies an anisotropic fluid, requires the introduction of a 
surface tension, for allowing the application of the algorithm. 
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C. Tolman VI type model 



Our last example is based on the Tolman VI solution. Accordingly the effective variables 
for this model will be 



r 

and 



^ (100) 



p _ g(9- bK{r/r^)) 

(9-6(r/r E ))r 2 ' K ' 

where g and b are functions of a, to be obtained from fl9~T|). Then, 

3(1 — F) , , 

P=- o ■ 102 

y 24vrr 2 v ; 

Using fl?4] ) and ( |75|) we get 

m = m 2 r/r s , (103) 



«/ = ZnF + S -f {4/n(r/r E ) + 8/n (^^) } " ^ 

Finally, solving the surface equations for this model, m and v are completely determined 
and all physical variables can thereby be calculated. Besides the intrinsic physical interest 
of the equation of state of this "seed" model mentioned before, it is interesting because of 
the fact that the static limit of the model (unlike the previous ones) is "unstable", in the 
sense that it requires a specific value of the gravitational potential at the boundary, namely 
^s(0)/rs = 3/14. For values above (below) this, the sphere starts to collapse (expand). 

Figures (7) and (8) display the evolution of velocity (u) for different regions of the 
sphere, and for initial values of F corresponding to values of m^(0)/r^, above and below the 
equilibrium value, respectively. Figures (9) and (10), represent the evolution of the Weyl 
tensor (W), for some internal region and the boundary surface respectively, and initials 
values of F corresponding to values of ms(0)/rs, above and below equilibrium. Finally, 
figures (11) and (12) exhibit the behaviour of the shear (a) for different regions and initial 
values of F corresponding to values of ms(0)/rs, above and below equilibrium. 

We shall comment on these graphics in the next Section. 



V. CONCLUSIONS 

A method has been presented, which allows for the description of radiating selfgravitating 
relativistic spheres. In its most general form, the approach incorporates the two limiting 
cases of radiation transport (free streaming and diffusion) as well as the possibility of dealing 
with anisotropic fluids. 
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The cornerstone of the algorithm is an ansatz, based on an specific definition of the 
post-quasistatic approximation, namely: considering different degrees of departure from 
equilibrium, the post-quasistatic regime (i.e. the next step after the quasistatic situation) 
is defined as that, characterized by metric functions whose radial dependence is the same as 
that of the quasistatic regime. This in turn implies, that the effective variables defined above, 
share the same radial dependence as the correspondig physical variables of the quasistatic 
regime. The rationale behind this definition seems intelligible when it is remembered that in 
the latter case (the quasistatic) the effective variables share the same radial dependence as 
that of the physical variables in the static regime. Thus, starting with a static configuration, 
the first "level" off equilibrium, beyond the quasitatic situation, is represented by the post- 
quasistatic regime. 

Once the static ("seed") solution has been selected, then the definition of the effective 
variables together with surface equations, allows for determination of metric functions, which 
in turn lead to the full description of physical variables as functions of the timelike coordinate, 
for any region of the sphere. In this process, depending on the kind of matter and/or the 
prevailing transport approximation, additional equations of state and/or transport equations 
and/or some of the surface variables (e.g. the luminosity) have to be specified. 

All physical variables having been found (particularly the energy density and the radial 
pressure) then we may, in principle, go to the next step assuming that the effective variables 
now share the same radial dependence as that of the physical variables just obtained. In this 
sense the algorithm may be regarded as an iterative approach. For obvious reasons we have 
restrained ourselves to the first step of the process. It remains to be seen if available physical 
evidence justifies to go through the complexities associated with the "post-post-quasitatic" 
approximation. 

In order to illustrate the method, and without the pretension of modeling specific astro- 
physical scenarios, we have presented three examples, in the simplest (adiabatic) case. 

In the first model, the profiles of the shear and the Weyl tensor, clearly illustrate the "dy- 
namics" of the model, tending to zero in the static limit. The fact that these two quantities 
vanish in the quasitatic regime (for this specific model), brings out further their relevance in 
the treatment of situations off equilibrium. On the other hand however, the velocity profiles 
show almost no difference between the two regimes. Deviations from homology contraction 
due to relativistic gravitational effect are also indicated. 

The purpose of the second example was to illustrate the implementation of the algorithm, 
for anisotropic fluids. The very particular form of the "seed" equation of state of this model, 
imposses discontinuity (surface tension) of radial pressure at the boundary. Of course such 
discontinuity vanishes at the static (or quasistatic) regime. 

Finally, a model based on the Tolman VI solution was presented. This static solution, as 
was already mentioned, requires a specific value of ms(0)/rs, accordingly any deviation from 
this value leads to deviations from the static regime (observe that the quasitatic regime is 
incompatible with this solution). The velocity profiles indicate that all regions either expand 
or contract, and therefore, cracking (different signs of the velocity for different regions of the 
sphere) will not occur [j36 ]. This is consistent with the established fact that cracking only 



occurs for anisotropic fluids or isotropic fluids with outgoing radiation in the free streaming 
approximation. 

Also, the profiles of the Weyl tensor and the shear, clearly diverging from the initial 
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values as time proceeds and the evolution becomes more and more "dynamic" , stress once 
again their roles in describing departures off equilibrium. 
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FIG. 1. Energy denstity, pm(0) 2 , as a function of (dimenssionless) time (a) for the 
Schwarzschild-type model. The initial conditions are: A(0) = 5, F(0) = 0.6 and $1(0) = —0.1. 
Curves represent different regions: r/rs = 0.25 (continuous line); 0.50 (dashed line); 0.75 
(short-dashed line) and 1.00 (dotted line). 
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FIG. 2. Radial pressure, P r m(0) 2 , as a function of time for the Schwarzschild-type model. The 
initial conditions are: A(0) = 5, F(0) = 0.6 and f2(0) = —0.1. Curves represent different regions: 
r/r-z = 0.25 (continuous line); 0.50 (dashed line); 0.75 (short-dashed line) and 1.00 (dotted line). 
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FIG. 3. Radial velocity, to, as a function of time for the Schwarzschild-type model. The initial 
conditions are: A(0) = 5, -F(O) = 0.6 and f2(0) = —0.1. Curves represent different regions: 
r/rs = 0.25 (continuous line); 0.50 (dashed line); 0.75 (short-dashed line) and 1.00 (dotted line). 



26 



0.03 



0.025 - 



0.02 - 



to 

CD 



CO 



0.015 - 



0.01 - 



0.005 1 1 1 ' 1 1 

2 4 6 8 10 

Time 

FIG. 4. Shear, am(0), as a function of time for the Schwarzschild-type model. The initial 
conditions are: ^4(0) = 5, -F(O) = 0.6 and O(0) = —0.1. Curves represent different regions: 
r/rs = 0.25 (continuous line); 0.50 (dashed line); 0.75 (short-dashed line) and 1.00 (dotted line). 
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FIG. 5. Weyl tensor, W/m(0), as a function of time for the Schwarzschild-type model. The 
initial conditions are: A(0) = 5, F(0) = 0.6 and f2(0) = —0.1. Curves represent different regions: 
r/rs = 0.25 (continuous line); 0.50 (dashed line); 0.75 (short-dashed line) and 1.00 (dotted line). 
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FIG. 6. Radial velocity, u, as a function of r/rs for a = 10. The initial velocity at the surface 
is —0.001. Curves represent different values of -F(O): 0.6 (continuous line); 0.96 (dashed line) and 
0.996 (short-dashed line). 
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FIG. 7. Radial velocity, uj, as a function of time for the Tolman Vl-type model. The initial 
conditions are: F(0) = 0.581428528 and O(0) = —0.0001. Curves represent different regions: 
r/rs = 0.2 (continuous line); 0.4 (dashed line); 0.6 (short-dashed line); 0.8 (dotted line) and 1.0 
(dot-dashed line). 
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FIG. 8. Radial velocity, uj, as a function of time for the Tolman-type model. The initial 
conditions are: F(0) = 0.561428547 and O(0) = —0.0001. Curves represent different regions: 
r/rs = 0.2 (continuous line); 0.4 (dashed line); 0.6 (short-dashed line); 0.8 (dotted line) and 1.0 
(dot-dashed line). 
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FIG. 9. Weyl tensor, W/m(0), at r/r^ = 0.4 as a function of time for the Tolman Vl-type 
model. The initial conditions are: F(0) = 0.581428528 (dashed line); F(0) = 0.561428547 (contin- 
uous line) and O(0) = -0.0001. 
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FIG. 10. Weyl tensor, W/m(0), at r/Vs = 1.0 as a function of time for the Tolman Vl-type 
model. The initial conditions are: F(0) = 0.581428528 (dashed line); F(0) = 0.561428547 (contin- 
uous line) and O(0) = -0.0001. 
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FIG. 11. Shear, am(0), as a function of time for the Tolman Vl-type model. The initial 
conditions are: -F(O) = 0.581428528 and O(0) = —0.0001. Curves represent different regions: 
r/rs = 0.2 (continuous line); 0.4 (dashed line); 0.6 (short-dashed line); 0.8 (dotted line) and 1.0 
(dot-dashed line). 
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FIG. 12. Shear, am(0), as a function of time for the Tolman-type model. The initial conditions 
are: -F(O) = 0.561428547 and O(0) = —0.0001. Curves represent different regions: r/r^ = 0.2 
(continuous line); 0.4 (dashed line); 0.6 (short-dashed line); 0.8 (dotted line) and 1.0 (dot-dashed 
line) . 
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